Microkinetic Analysis of the Oxygen Evolution Performance at Different Stages of Iridium Oxide Degradation

The microkinetics of the electrocatalytic oxygen evolution reaction substantially determines the performance in proton-exchange membrane water electrolysis. State-of-the-art nanoparticulated rutile IrO2 electrocatalysts present an excellent trade-off between activity and stability due to the efficient formation of intermediate surface species. To reveal and analyze the interaction of individual surface processes, a detailed dynamic microkinetic model approach is established and validated using cyclic voltammetry. We show that the interaction of three different processes, which are the adsorption of water, one potential-driven deprotonation step, and the detachment of oxygen, limits the overall reaction turnover. During the reaction, the active IrO2 surface is covered mainly by *O, *OOH, and *OO adsorbed species with a share dependent on the applied potential and of 44, 28, and 20% at an overpotential of 350 mV, respectively. In contrast to state-of-the-art calculations of ideal catalyst surfaces, this novel model-based methodology allows for experimental identification of the microkinetics as well as thermodynamic energy values of real pristine and degraded nanoparticles. We show that the loss in electrocatalytic activity during degradation is correlated to an increase in the activation energy of deprotonation processes, whereas reaction energies were marginally affected. As the effect of electrolyte-related parameters does not cause such a decrease, the model-based analysis demonstrates that material changes trigger the performance loss. These insights into the degradation of IrO2 and its effect on the surface processes provide the basis for a deeper understanding of degrading active sites for the optimization of the oxygen evolution performance.


Model Implementation
The microkinetic model is implemented and simulated in MATLAB 2020b. The set of differential equations is solved by applying an ode23s solver. High accuracy of the results is ensured by setting the absolute error tolerance to 10 -12 and the maximum step size to 0.01 s.

Model Parameters
Parameter identification is crucial in order to gain a valid microkinetic model and, thus, valid model output, such as physically meaningful merits and insights into the dynamic behavior of the simulated system. For our study one can divide the parameters into three main categories: 1) Parameters of which the values are defined by the model assumptions e.g. the properties gained directly from the presumed mechanism itself. 2) Experimentally accessible parameter values, which can be measured easily or quantified by analyzing experimental data e.g. temperature and electrolyte resistance. 3) Parameters which are neither accessible directly with experimental methods nor defined by the model assumptions e.g. density of active sites or thermodynamic energies. These demand for either ab initio/DFT calculations or model-based identification algorithms to quantify their value from experiments. Unfortunately, DFT calculations show strong discrepancies in the resulting reaction energy values depending the computational details and assumptions 1 , and they usually have high uncertainties, which in extreme cases may amount up to 0.66 eV. In addition, mostly reaction free energies were published so far, but few activation barrier energies, and to the best of our knowledge also no interaction energies of the adsorbed species.
Hence, model-based parameter identification algorithms are required in order to completely parametrize the microkinetic model. In the following paragraphs, we will introduce the parameter identification process for all model parameters, which are needed to gain a valid model for dynamic simulations and further analyses.
The parameters of the first category are straightforward derived from the presumed mechanism.
The parameters of the second category are quantified by the experimental properties and results.
All experiments were conducted at a room temperature of T = 25 °C, and the geometrical electrode area was A = 0.1963 cm 2 . The electrolyte resistance of R = 18 Ω was determined by electrochemical impedance spectroscopy measurements; the value was set to the impedance obtained at a phase angle of φ = 0° in the high frequency range between f = 10 5 Hz and 10 4 Hz.
The activity of electrolyte species at the electrode interface was assumed to equal the respective bulk value due to fast transport compared to the consumption or production rates by the reactions.
This assumption is reasonable, since in cyclic voltammetry experiments an increasing thickness of the Nernst diffusion layer accomplished by applying slower rotation speed indicated no significant limitation. Therefore, the activities were implemented as constant values of ! = 1 for water and " = 0.2 for protons, as the 0.1 M sulfuric acid is known to dissociate completely at low concentrations. algorithm, which is described to full extend in the following subsection. Besides the density of active sites ρ and the double layer capacitance Cdl, this procedure was also applied to estimate the reaction free energy Δ& ' and activation free energy Δ& ( values of all seven reactions and the interaction free energy values of all seven adsorbed species Δ& )*+ .
It is noteworthy to state here, that in contrast to both the transition state theory, which is derived for forward reactions only, and the more generalised Marcus theory 2 , we do not assume a preexponential frequency factor other than k0 = 343.2 s -1 . This simplification is valid since we model the electrocatalytic system at constant temperature T and, hence, the exponential relation between the frequency factor k0 and the activation energy ΔGa makes it unfeasible to identify both parameters independently. For further explanation see the following section.

Parameter Identification Process
Figure S1: Schematic illustration of the two-step parameter identification process including the global parameter search of 10 6 randomly selected sets of parameters and the local parameter optimization of the best 250 sets.
To estimate the parameters of the third category, which are the reaction free energies, the activation free energies and the interaction free energies, the double layer capacitance and the density of active sites, we conducted an optimization procedure in two steps, as illustrated in Figure S1 a).
In a first step, the dynamic simulations were performed by randomly defining 10 6 sets of parameters within suitable ranges shown in Table S1. The ranges of the reaction free energy values of the electrochemical steps were set in such way that they were able to reproduce the features in the CV e.g. peak positions. For chemical steps, the ranges of the reaction free energies were selected in between wider limits, as they are not easily specifiable by the experimental current response. The reaction free energy of the oxygen evolution (step 7: *OO → * + O2) is constrained by the thermodynamic formation energy ∑Δ& ' -= 4.92 eV of the overall OER (2 H2O → 2 H2 + O2). The preexponential frequency factor k0 was set by the analysis of the peak-to-peak potential of anodic and cathodic redox transition, which converges to 0 V at a value of roughly 2 -= 343.2 s 5 when assuming no activation barrier as shown in Figure S2. This value was chosen because we identify the activation free energy from the peak-to-peak potential. In consequence, the lower limits of the activation free energies were defined as Δ& ( = 0 eV. We assume as the upper range a value of 0.15 eV, since no significant peak to peak difference is observable in the CVs of the pristine material. One exception is the oxygen detachment (step 7: *OO → * + O2) which was found by DFT calculations to have an activation energy up to 0.58 eV 3,4 and, thus, the respective activation free energy upper limit was increased to 1.65 eV. Lower and upper limits of the interaction energies of adsorbed species were set to 0 eV and 0.2 eV, respectively. The density of active sites and the double layer capacitance were found suitable to describe the experimental data in the ranges of [8·10 -5 , 20·10 -5 ] mol m -2 and [18, 20] F m -2 respectively.   Table 1, respectively.
The pattern search algorithm was also applied to identify the model parameter values of the steadily degrading catalytic system. Successively, the CV curves which were measured after each 30 minutes of operation at 1.6 V were used as input data =>? in the objective function to minimize the rmse. In every step, the values from the previously optimized set served as initial parameters.
The simulation results as well as the rmse values are shown in the main manuscript in Figure 7.

Electrochemical Degradation Results
Application of a constant operating potential for a total duration of 15 hours in the acidic electrolyte leads to a decrease in the OER activity of the catalytic system. The dependence of changes in the CV on the constant potential is shown in Figure S3  An additional experiment with constant potential of 1.6 V applied for 15 hours followed by a constant potential of 1.2 hours for 5 hours was conducted. The OER activity was measured by CV and is shown in Figure S4. After the performance decrease during operation at 1.6 V, a slight recovery of the activity is visible while applying a potential of 1.2 V afterwards.

CV Simulations with Different Scan Rates
Dynamic microkinetic simulations allow for arbitrary change of the model input, i.e. the applied potential over time. With the aim to prove that our model is able to reproduce the typical experimental variation in the potential scan rate, we show the simulated and the measured CV curves in Figure S5 a)

Profile root mean square error analysis
To ensure parameter identifiability of the reaction energy values, profile root mean square error (rmse) analysis was conducted. The method follows the profile likelihood analysis described elsewhere, 5 Table S2.
To simulate the overall OER, a high amount of 23 parameters is needed to be identified. To parameter defines an identifiable feature in the dynamic current response.  Bandarenkas publication, a simple equivalent circuit model, was used to reproduce the spectra and identify parameters. 6 Thus, the adsorption capacitance Ca was quantified as given in Figure S8 b).
It is the only parameter correlated to the ECSA. It shows the same increasing trend with increasing potential up to 1.6 V as reported for IrOx. 6 To determine the ECSA from Ca, we used the procedure and value of the active-area specific adsorption capacitance Ca' from the Bandarenka group 6 . Here, the capacitance Ca,B measured by EIS was normalized to the ECSAB of an IrOx thin film determined by atomic force microscopy, which resulted in Ca' = Ca, B / ECSAB = 135 ± 25 µF cm -2 . This value is used to calculate the actual ECSA of the present system ECSA = Ca / Ca' = 2.7 ± 0.6 cm 2 .

Oxidation State Validation
The in-situ X-ray cell built for this analysis is based on the development of Binninger et al. 7 The cell preparation is explained to full extend in one of our previous reports from Czioska et al. 8  The raw data was then treated using Athena and Artemis. 9 The energies of the spectra were calibrated and aligned using reference metal foil channels to correct any energy shifts, and the spectra were normalized. The resulting position of the white line is shown in Figure S9 and taken as an indicator for the change in Ir oxidation state. Although XAS radiation probes not only the active surface material but also the bulk material, a significant shift is observable with potential.
Since the bulk material does not change significantly, the shift can be attributed to the Ir atoms on the outer surface. With the aim to compare the results to the model output we explain how we calculate the mean oxidation state in the following paragraph.
The surface coverage, which corresponds to the share of surface covered by a species or its probability to be found at the surface, is multiplied with the difference in oxidation state. The results are shown in Figure S9 for direct comparison to the operando XAS measurements.
In direct comparison to the white line peak position of the Ir L3-edge absorption spectra at constant applied potentials, simulated and experimental results correspond on two main trends: 1) For both, an increase in the mean oxidation state with potential is observed in the range of 0.6 V to 1.35 V.
2) Both also show that at higher potentials, a maximum oxidation state is reached, after which the oxidation state decreases, thus, the active sites are reduced. While optimizing the parameters during the steady degradation of the catalytic system, the root mean squared error (rmse) is tracked and given in Figure S11. The value of the density of active sites is always optimized with parameter combinations tested as denoted in the bottom right of each graph to account for material dissolution, particle detachment and loss of binder material.
The change in proton concentration, electrolyte resistance and double layer capacitance are individually evaluated and the maximum rmse values of 0.33 mA cm -2 , 0.18 mA cm -2 and 0.26 mA cm -2 are gained respectively in Figure S11 a-c). Although the change in electrolyte resistance is evaluated with lowest rmse so far, a maximum parameter value of R = 392 Ω is required to describe the completely degraded system, which is in contrast to the experimentally observed 18.2 Ω. Due to the high discrepancies we can disprove any correlation between the electrolyte-related parameters and the degradation behavior.
The free energy parameters are analyzed with the same method, and the results are shown in Figure S11 d-i). Low values of the rmse and, thus, good reproduction of the experimental behavior is gained by optimizing all activation free energy values in Figure S11 d). We found the activation free energy to be the best descriptor for the degradation process with a maximum rmse value not exceeding 0.11 mA cm -2 , which is significantly lower compared to the error values of 0.18 mA cm -2 for the reaction energy in Figure S11 e) and 0.22 mA cm -2 for the interaction free energies in Figure S11 f), respectively. Optimizing only activation free energies of the deprotonation steps hardly increase the rmse in Figure S11 g). Also, the change of the reaction energy in Figure S11 h) or the interaction energy in Figure S11 i) of only deprotonation steps of gives comparable rmse values with a change in all steps. Evaluated with a high rmse value of up to 0.55 mA cm -2 in Figure   S11 j), the change of the density of active sites is not able to explain the experimentally observed degradation behavior solely. Best reproduction is obtained by the optimization of the activation and reaction free energies of the deprotonation steps, evaluated with a rmse remaining for each evaluated timestep below the 0.1 mA cm -2 threshold in Figure S11 k). The simulation results and identified parameter values are shown in the main document in Figure 7.